Update notes v12

Update notes v11

Load libraries and define parameters

NB: this Rmd can be used for all types of ratings: emotion, arousal, valence

Abbreviations: - sub : subjects(s) - rat : rating(s)

Load libraries

library(tidyverse)
library(future)
library(furrr)
library(tictoc)
library(RNifti)
library(proxy) # distances
library(profvis)
library(DT)
library(formattable)
library(janitor)
library(ComplexHeatmap)
library(circlize)  # For color mapping functions
library(grid)  # For gpar()
# library(pheatmap)
# library(heatmaply)
library(papayaWidget)
library(ggstatsplot)


source("funs_V12_ROIs.R")

Copes/emotion dictionary specific for emotion_prediction

metadata_copes <- tribble(
  ~cope, ~emotion,
  1, "anger",
  2, "disgust",
  3, "fear",
  4, "happy",
  5, "neutral",
  6, "pain",
  7, "sad"
)

Options for different flavours

# Do you want to remove neutral cope(s) from the RSA analysis?
# If "YES", the results filename will have option _RM (remove neutral)
remove_neutrals = "NO"

# choose the type of betas to be used following this correspondence table:
# rsa_ROI : one_ev_per_movie[_minus_neutral]
# rsa_emotion_predictors_ROI : emotion_predictors[_minus_neutral]
# rsa_emotion_high_low_predictors_ROI : emotion_high_low_predictors[_minus_neutral]
copes_type <- "emotion_predictors_minus_neutral"

# Choose one set of ROIs - e.g. an atlas - from  ${bd_atlases} (ROIS_REPO)
atlas_filename <- "Yeo17.nii.gz"

# Schaefer.nii.gz
# Schaefer.nii.gz
# Yeo17.nii.gz
# Yeo17.nii.gz

# Atlases repos:
# https://www.lead-dbs.org/helpsupport/knowledge-base/atlasesresources/cortical-atlas-parcellations-mni-space/
# https://github.com/neurodata/neuroparc
# https://www.fmrib.ox.ac.uk/datasets/brainmap+rsns/

Define paths and rsa_[flavour] specific option

bd="/data00/leonardo/RSA/analyses"
bd_ratings = paste0(bd,"/RATINGS")
bd_atlases = paste0(bd,"/ROIS_REPO")
bd_results = paste0(bd,"/RSA_ROI_APP/results_RSA_ROI")


ratings_type <- c("emotion","arousal","valence","aroval")


rsa_flavour="rsa_emotion_predictors_ROI"


# Choose the distance metric to use for fmri RDM
# The metric for ratings RDM *should be* euclidean, since arousal and valence
# ratings have only one value
# Supported methods: 'pearson', 'spearman', 'euclidean', 'cosine' or 'mahalanobis'
dist_method_rating <- "euclidean"
dist_method_fmri = "euclidean"
dist_method_rsa = "pearson"

# Vector of zeropadded sub_ids
subs_file <- "/data00/leonardo/RSA/sub_list.txt"
subs <- sprintf("%02d",readLines(subs_file) %>% as.numeric)

# ---------- ONLY TOP RATERS BELOW ------------
subs <- c("02","03","12","11","22","09","29","28","26","32","23","15","20","19")

Create a df_path_copes with the location of the 56 cope niis from the one_movie_per_ev model

Extract the pathname of all copes using the list.files() function. Also define a copes_numba vector with all the copes numbers.

NB: The cope numbers in the cope column are NOT zeropadded since this is how they come out from FSL Feat

# load df_path_copes and add metadata (label and emotion/neutral)
df_path_copes <- import_df_path_copes(bd, copes_type, rsa_flavour) %>% 
  inner_join(metadata_copes, by="cope") %>% 
  relocate(sub, cope, emotion, path)

# If remove_neutrals is "YES", remove neutral copes and arrange by sub/cope (just to be sure) 
if (remove_neutrals == "YES") {
  df_path_copes <- df_path_copes %>%  
  filter(emotion != "neutral") %>% 
  arrange(sub,cope)
} 

# ---------- SPECIFIC TO EMOTION PREDICTION -------------
# add a label colum to the df_path_copes which represents the emotion
# (for the original 56 movies it was high_low_code, e.g. JK_A_high)

df_path_copes$label <- df_path_copes$emotion

# ---------- SPECIFIC TO EMOTION PREDICTION -------------



copes_numba <- df_path_copes$cope %>% unique

cat(
  "there are",length(copes_numba)," copes including ",
  nrow(df_path_copes %>% filter(sub==subs[1], emotion=="neutral"))," neutral copes \n"
)
## there are 7  copes including  1  neutral copes
# df_path_copes
# length(copes_numba)

df_path_copes
## # A tibble: 182 × 5
##    sub    cope emotion path                                                label
##    <chr> <dbl> <chr>   <chr>                                               <chr>
##  1 02        1 anger   /data00/leonardo/RSA/analyses/emotion_predictors_m… anger
##  2 02        2 disgust /data00/leonardo/RSA/analyses/emotion_predictors_m… disg…
##  3 02        3 fear    /data00/leonardo/RSA/analyses/emotion_predictors_m… fear 
##  4 02        4 happy   /data00/leonardo/RSA/analyses/emotion_predictors_m… happy
##  5 02        5 neutral /data00/leonardo/RSA/analyses/emotion_predictors_m… neut…
##  6 02        6 pain    /data00/leonardo/RSA/analyses/emotion_predictors_m… pain 
##  7 02        7 sad     /data00/leonardo/RSA/analyses/emotion_predictors_m… sad  
##  8 03        1 anger   /data00/leonardo/RSA/analyses/emotion_predictors_m… anger
##  9 03        2 disgust /data00/leonardo/RSA/analyses/emotion_predictors_m… disg…
## 10 03        3 fear    /data00/leonardo/RSA/analyses/emotion_predictors_m… fear 
## # ℹ 172 more rows

Read atlas and create a vector of atlas labels

atlas_path <- paste0(bd_atlases,"/",atlas_filename)

atlas_nii <- readNifti(paste0(bd_atlases,"/",atlas_filename))
region_labels <- atlas_nii[atlas_nii > 0] %>% unique %>% sort

Calculate RATINGS RDMs

# NB: ratings_type is now a *vector*, e.g.: 
# ratings_type <- c("emotion","arousal","valence","aroval")
# to test the fn below for one rating type : do_RDMs_rats("emotion")


do_RDMs_rats <- function(ratings_type, remove_neutrals) {

  ratings_path <- paste0(bd_ratings,"/",ratings_type,"_ratings.csv")
  rats <- read_csv(ratings_path)
  
  
  # ------- IMPORTANT : AVG ACROSS MOVIES FOR EMOTION PREDICTORS MODEL -----
  rats <- rats %>% 
    select(sub, emotion, starts_with("r_")) %>%
    group_by(sub, emotion) %>% 
    reframe(across(starts_with("r_"), median, na.rm = TRUE)) %>% 
    ungroup %>% 
    arrange(sub, emotion)
  
  # ------- IMPORTANT : AVG ACROSS MOVIES FOR EMOTION PREDICTORS MODEL -----
  
  
  if (remove_neutrals == "YES") {
    rats <- rats %>% filter(emotion != "neutral")
  }
  
  rdm <- rats %>%
    filter(sub %in% subs) %>%
    select(sub, starts_with("r_")) %>%
    group_by(sub) %>%
    nest %>%
    mutate(!!paste0("RDM_",ratings_type) := data %>% map(~ DDOS_vec(.x, dist_method_rating))) %>%
    ungroup %>% 
    # remove the sub with the prospect of merging rdms of different ratings/models
    select(!data)
  
  return(rdm)
  
}

# calculate rdms for all ratings_type and return one column for each ratings_type
RDMs_rats <- map_dfc(ratings_type, ~ do_RDMs_rats(.x, remove_neutrals)) %>%
  # remove the duplicated sub_ columns (generated by do_RDMs_rats) 
  janitor::clean_names() %>% 
  mutate(sub = sub_1) %>% 
  select(sub, starts_with("rdm"))
## Warning: There was 1 warning in `reframe()`.
## ℹ In argument: `across(starts_with("r_"), median, na.rm = TRUE)`.
## ℹ In group 1: `sub = "02"` and `emotion = "anger"`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))
# # example heatmap
# plot_tril(
#   RDMs_rats[1,]$rdm_emotion, model = "emotion",
#   fontsize=10, side_mm=80, reord = FALSE
# )

Calculate FMRI RDMs

IMPORTANT : Since we already removed neutral from df_path_copes we don’t need to filter anything anymore here.

# Fn to calculate the rdm_fmri for one sub
# Test for one sub: do_RDM_fmri(subs[1], copes_numba, df_path_copes)

do_RDM_fmri <- function(sub_id, copes_numba, df_path_copes) {
  
  # load all the copes nii.gz into a df
  df_copes <- load_sub_copes(sub_id, copes_numba, df_path_copes)
  
  # the code below does the calculation of the tril of RDM_fmri for all rois for one sub
  
  one_sub_RDM_fmri <- tibble(
    sub = sub_id,
    roi = region_labels
  ) %>%
  
    # extract idx_roi for atlas voxels in that roi
    mutate(idx_roi = roi %>% map(~ which(atlas_nii == .x)) ) %>%
    
    # extract the df_copes for the voxels in that region (idx_roi)
    mutate(df_copes_region = idx_roi %>% map(~ df_copes[.x,]) ) %>%
    
    # calculate tril of RDM_roi (output as numeric vector)
    mutate(rdm_fmri = df_copes_region %>% map(~ DDOS_vec(t(.x), dist_method_fmri)) ) %>% 
    
    # select only sub, roi (numba) and RDM_roi
    select(sub, roi, rdm_fmri) 
  
  return(one_sub_RDM_fmri)
}



# Calculate RDMs_fmri for all subs

plan(multisession, workers = 5)

RDMs_fmri <- subs %>% future_map_dfr(
  ~ {
    paste0("Calculating RDMs for sub ",.x,"\n") %>% cat
    do_RDM_fmri(.x, copes_numba, df_path_copes)
  }
)
## Calculating RDMs for sub 02
## Calculating RDMs for sub 03
## Calculating RDMs for sub 12
## Calculating RDMs for sub 11
## Calculating RDMs for sub 22
## Calculating RDMs for sub 09
## Calculating RDMs for sub 29
## Calculating RDMs for sub 28
## Calculating RDMs for sub 26
## Calculating RDMs for sub 32
## Calculating RDMs for sub 23
## Calculating RDMs for sub 15
## Calculating RDMs for sub 20
## Calculating RDMs for sub 19
plan(sequential)

# RDMs_fmri

Calculate RSA

# The dist_method_rsa is set above.
# can be pearson, spearman
# dist_method_rsa <- "spearman"

# Join fmri and rats RDMs and create a copy of the rat RDM for each sub/roi
RSA <- right_join(RDMs_fmri, RDMs_rats, by = "sub") %>%
  group_by(sub) %>%
  
  # Dynamically create rsa_ columns for every ratings type
  mutate(across(
    all_of(paste0("rdm_", ratings_type)),
    ~ map2_dbl(.x, rdm_fmri, ~ cor(.x, .y, method = dist_method_rsa)),
    .names = "rsa_{.col}"
  )) %>%
  mutate(across(starts_with("rsa_"), round, 3)) %>%
  ungroup()
## Warning: There were 17 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `across(...)`.
## ℹ In group 6: `sub = "15"`.
## Caused by warning in `cor()`:
## ! the standard deviation is zero
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 16 remaining warnings.
# Get the mean RSA for each ratings_type
RSA_mean <- RSA %>%
  select(roi, starts_with("rsa_")) %>%
  group_by(roi) %>%
  reframe(
    across(
      starts_with("rsa_"),
      list(mean = ~ round(mean(.x, na.rm = TRUE), 2))
    )
  ) %>% 
  rename_with(
    ~ str_replace(.x, "rsa_rdm_", ""),
    starts_with("rsa_")
  )

RSA_mean %>% datatable()

Calculate mean tril for ratings and for fmri (for each roi in the fmri case)

# Takes a df and a "column" containing a nested vector in each row
# and returns the element-wise mean of that vector across rows
get_elementwise_mean <- function(mydf, nested_col, fun=mean) {
  m <- matrix(
    mydf[[nested_col]] %>% unlist,
    nrow = nrow(mydf),
    ncol = length(mydf[[nested_col]][[1]]),
    byrow = TRUE
  )
  
  m_summary <- apply(m, 2, fun)
  return(m_summary)
}


# mean RDMs_rats
mean_RDMs_rats <- RDMs_rats %>% select(starts_with("rdm_")) %>% colnames %>%  
  map_dfc(~{
    mean_tril <- get_elementwise_mean(RDMs_rats, .x, fun = mean)
    tibble(!!.x := mean_tril)
  })



# mean RDMs_fmri
plan(multisession, workers = 5)

mean_RDMs_fmri <- RDMs_fmri$roi %>% 
  future_map_dfr(~{
    mean_tril <- get_elementwise_mean(RDMs_fmri %>% filter(roi == .x), "rdm_fmri")
    tibble(roi = .x, mean_tril = mean_tril)
  }) %>% 
  group_by(roi) %>% 
  nest

plan(sequential)

Plot mean RDMs_rats

# Plot mean RDMs_rats
ratings_type %>% 
  map(
    ~ plot_tril(
      mean_RDMs_rats[[paste0("rdm_",.x)]], model = .x,
      fontsize=7, side_mm=80, reord = FALSE
    )
  )
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Plot mean RDMs_fmri for all ROIs

mean_RDMs_fmri
## # A tibble: 17 × 2
## # Groups:   roi [17]
##      roi data              
##    <dbl> <list>            
##  1     1 <tibble [294 × 1]>
##  2     2 <tibble [294 × 1]>
##  3     3 <tibble [294 × 1]>
##  4     4 <tibble [294 × 1]>
##  5     5 <tibble [294 × 1]>
##  6     6 <tibble [294 × 1]>
##  7     7 <tibble [294 × 1]>
##  8     8 <tibble [294 × 1]>
##  9     9 <tibble [294 × 1]>
## 10    10 <tibble [294 × 1]>
## 11    11 <tibble [294 × 1]>
## 12    12 <tibble [294 × 1]>
## 13    13 <tibble [294 × 1]>
## 14    14 <tibble [294 × 1]>
## 15    15 <tibble [294 × 1]>
## 16    16 <tibble [294 × 1]>
## 17    17 <tibble [294 × 1]>
seq(nrow(mean_RDMs_fmri)) %>% 
  map(
    ~ plot_tril(
      mean_RDMs_fmri[.x,]$data, model = paste0("fmri_roi_",.x),
      fontsize=7, side_mm=80, reord = FALSE
    )
  )
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

## 
## [[17]]

Elements for the app

RSA table

RSA_mean %>% datatable()

Boxplot for one roi

roi_numba = 1

RSA %>% 
  select(sub, roi, starts_with("rsa_")) %>% 
  filter(roi == roi_numba) %>% 
  select(starts_with("rsa_")) %>% 
  pivot_longer(cols = starts_with("rsa_"), names_to = "model") %>%
  mutate(model = str_replace(model, "rsa_rdm_","")) %>% 
  ggwithinstats(
    x = model, y = value,
    type = "p",   # p, np, r, b
    results.subtitle = TRUE
  )

View one ROI

prepare_papaya <- function(atlas_filename, roi_numba) {
  
  f <- function(nii_filename) {
    return(paste0(bd_atlases,"/", nii_filename))
  }
  
  nii_atlas <- RNifti::readNifti( f(atlas_filename) )
  # view(nii_atlas)
  
  nii_atlas_thr <- ifelse(nii_atlas == roi_numba, 1, 0)
  RNifti::writeNifti(nii_atlas_thr, template = nii_atlas, paste0(bd_atlases,"/tmp_ROI.nii.gz"))
  
  papaya(
    c(f("Dummy.nii.gz"), f("MNI152_T1_2mm_brain.nii.gz"), f("tmp_ROI.nii.gz")),
    options = list(
      papayaOptions(alpha = 1, lut = "Grayscale"),
      papayaOptions(alpha = 0.5, lut = "Grayscale"),
      papayaOptions(alpha = 0.5, lut = "Red Overlay", min = 0, max = 5)
    ),
    interpolation = FALSE,
    orthogonal = FALSE
  )
  
}

prepare_papaya(atlas_filename, roi_numba = 6)

Plot MEAN tril for fmri and all models for one roi

global_reord = FALSE
global_fontsize = 7
global_side_mm = 100

roi_numba = 6

# fmri one roi
tril_fmri <- plot_tril(
  mean_RDMs_fmri[roi_numba,]$data, model = paste0("fmri_roi",roi_numba),
  fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)


# emotion
tril_emotion <- plot_tril(
  mean_RDMs_rats$rdm_emotion, model = "emotion",
  fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)

# arousal
tril_arousal <- plot_tril(
  mean_RDMs_rats$rdm_arousal, model = "arousal",
  fontsize = global_fontsize, side_mm = global_side_mm, reord = global_reord
)

tril_fmri + tril_emotion

Save results in RSA_ROI_APP

# We save the following into /data00/leonardo/RSA/analyses/RSA_ROI_APP
# RSA
# RSA_mean
# atlas_filename
# mean_RDMs_fmri
# mean_RDMs_rats


create_filename_results <- function(atlas_filename, copes_type, remove_neutrals) {
  atlas_root <- str_replace(atlas_filename, ".nii.gz", "")

  if(str_detect(copes_type,"one_ev_per_movie")) {EV_abbr <- "_OEPM"}
  if(str_detect(copes_type,"emotion_predictors")) {EV_abbr <- "_EP"}
  if(str_detect(copes_type,"emotion_high_low_predictors")) {EV_abbr <- "_EHLP"}
  
  option_MN <- ifelse(str_detect(copes_type,"_minus_neutral"),"_MN","")

  option_NN <- ifelse(remove_neutrals == "YES", "_RN", "")
  
  results_filename <- paste0(atlas_root, EV_abbr, option_MN, option_NN)
  return(paste0(bd_results,"/",results_filename,".RData"))
  
}

results_file <- create_filename_results(atlas_filename, copes_type, remove_neutrals)
results_file
## [1] "/data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI/Yeo17_EP_MN.RData"
# save the results
save(RSA, RSA_mean, atlas_filename, mean_RDMs_fmri, mean_RDMs_rats, df_path_copes, subs, file = results_file)

# load them again
# load(results_file)
# load("/data00/leonardo/RSA/analyses/RSA_ROI_APP/results_RSA_ROI/Yeo17_OEPM_MN_RN.RData")